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Abstract 

We study the dynamics of ultracold Bosons in a double-well potential within the two-mode 
Bose-Hubbard model by means of semiclassical methods. By applying a WKB quantization we 
■ find analytical results for the energy spectrum, which are in excellent agreement with numerical 

bJQ! 

exact results. They are valid in the energy range of plasma oscillations, both in the Rabi and 
' the Josephson regime. Adopting the reflection principle and the Poisson summation formula we 

derive an analytical expression for the dynamics of the population imbalance depending on the 
few relevant parameters of the system only. This allows us to discuss its characteristic dynamics, 
especially the oscillation frequency, and the collapse- and revival time, as a function of the model 
parameters, leading to a deeper understanding of Josephson physics. We find that our fomulae 
match previous experimental observations. 

PACS numbers: 03.65.Sq, 03.75.Lni, 05.30.-d 
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I. INTRODUCTION 



Fundamental issues of non-equilibrium physics of interacting many-body quantum sys- 
tems and of phase coherence and phase stability, in particular, have a long history. A simple 
yet relevant model, the two-site Bose-Hubbard Hamiltonian, features phase and fluctua- 
tion decay, and also revivals and thus, over the years, many thorough investigations of its 
quantum dynamics have appeared. Most remarkably, recent experiments involving ultracold 
Bose gases trapped in an effectively one-dimensional double-minimum potential represent 
an almost ideal realizations of this fundamental model [l, 2!, with the fascinating possibility 
to vary relevant model parameters over a wide range. 

A full many-body calculation of the dynamics of an interacting, trapped ultracold Bose 
gas is only possible for a very small number of particles, even for weakly interacting Bosons. 
Most often a mean field approximation in form of the Gross-Pitaevskii equation is applied, 
which provides good results for low temperatures and for a large number of particles N, if 
only for a limited time and a limited set of observables. These limits are intensively studied. 
Once the field operators are replaced by a c-number field, some truly quantum phenomena, 
e.g. wave function revivals, cannot be described. The double-well potential provides an 
ideal playground to analyze these issues. Thus, a purely classical field approach quickly 
comes to its limits, and the question arises whether semidassical methods can improve 
the theoretical treatment of such bosonic systems, allowing us in the future to study more 
challenging problems whose many-body Schrodinger equation can no longer be solved fully 
numerically. 

A number of articles deal with the discussion of the consequences of the mean-field ap- 
proximation and many-body quantum corrections Q] and the many-body quantum and 



classical dynamics in phase s 



the double-well system. In [6 
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^furthermore, semidassical methods were applied to 
B|-|25| a WKB quantization is adopted to analyze the 
energy spectrum and the wave functions in certain parameter regions. 

Despite this fair amount of investigations, it is remarkable to realize that - leaving some 
fairly straightforward cases aside - no analytical expressions for the relevant dynamical 
quantities appear to be known. Thus, the purpose of this article is to find a generally 
applicable analytical description of the population imbalance dynamics of an ultracold Bose 
gas in a double-well potential by applying semidassical methods. Since the full quantum 
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dynamics can be determined numerically up to many thousand particles, we are able to 
compare to exact results. Clearly, the interesting case of very large N ^ oo can no longer 
be investigated numerically, yet our analytical approach is suited to study this very limit in 
detail. 

At low temperatures a Bose-Einstein condensate in a double-well potential can be de- 
scribed by a two-mode approximation. The corresponding second quantized many particle 
two-site Bose-Hubbard Hamiltonian is written as 

Hbh = —T (^a\a2 + alai^ + U (^a[a[aiai + 0,202^^20.2 j + 5 {ni — 712) (1) 

with the creation and annihilation operators for a boson in the ith. well denoted by a|, 
hi with [oj, aj] = 5ij. Thus, the particle number operator of the iih. site is hi = a\ai. f/ is a 
measure for the on-site two-body interaction strength, T is a tunneling amplitude, which in 
the experiments can be controlled by varying the barrier hight. The tilt parameter 5 leads 
to an asymmetry in the one-particle site energies of the two wells and is used to initiate 
the dynamics. Note that in the standard notation adapted in Josephson physics we have 
Ej = NT and Ec = 4JJ 



It has been shown that the Bose-Hubbard Hamiltonian describes the dynamics of the 



bosons in the double- well potential properly js], provided that the interaction energy 
U is small compared to the level spacing of the trap potential, such that only the two 
lowest lying modes have to be taken into account. Transverse modes should also be 
suppressed. It should be mentioned that there are finer descriptions of the two-mode limit 
that also take into account tunnel coupling energies depending explicitly on the nonlinear 
two-body interaction term [9j]. In this work, however, we restrict ourselves to the standard 
Bose-Hubbard Hamiltonian (|T]). 



First, there are three qualitatively quite different regimes [10|, lUl with respect to 
crucial features of the energy spectra. They are best explained by introducing the 
parameter 

A ^ ^. (2) 

which thus separates the Rabi- (A < 1) from the so-called Josephson regime for which 
1 < A -C and the Fock regime with A ^ . 
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The Rabi regime is the non- interacting hmit A ^ 1, when the system consists of inde- 
pendent particles leading to an almost harmonic oscillator energy spectrum and thus, after 
an initial tilt, to plasma oscillations with the known plasma frequency u„ = 2Ty/l + A ^ 2T 

QQ. 

In the Fock regime all eigenenergies are grouped in doublets with a quasi-degenerate 
symmetric and antisymmetric state. Thus, the dynamics of the mean population imbalance 
follows an extremely slow evolution in time which is called self trapping. 

The Josepshon regime combines the two characteristics of the spectrum just discussed. 
We distinguish the self trapping regime E > 2NT from the plasma oscillating regime, where 
E < 2NT holds. In the former, the energy eigenstates appear as doublets again leading to 
self trapping. In the latter the energy eigenstates correspond to an (an-harmonic) oscillator 
spectrum and the population imbalance oscillates around zero. 

Thus, in the Josephson regime the dynamics will depend on the energy of the initial 
state. For low energies - the subject of this work - the dynamics undergoes plasma 
oscillations, for higher energies we see self-trapping, which is beyond the scope of this paper. 

In this article we have in mind an experiment as in reference so the double- well 
system is initially prepared in the ground state ipo of a tilted potential, i.e. 5 7^ in ([1]). 
Then, at t = it is quickly switched to a symmetric potential, i.e. 6 = 0. Starting from an 
initial population imbalance unequal to zero the system is left to evolve in time. 

In our paper we first discuss the spectrum using the semiclassical WKB- or Bohr- 
Sommerfeld quantization. We find a way to systematically obtain an approximate, useful 
expression for energies in the plasma oscillating regime. In order to describe imbalance dy- 
namics, we need to explore overlap matrix elements in the following section, which we do 
with the help of the reflection principle. We then apply the Poisson summation formula, 
which has a long history in semiclassical approaches to quantum dynamics. As a result, 
we find a useful expression for the time evolution of the imbalance, containing parameters 
that can be obtained analytically on the basis of the classical Hamiltonian. We then com- 
pare exact calculations with our new formula and find remarkable agreement over the whole 
relevant range of A, covering the known Rabi- but also the plasma oscillating Josephson re- 
gion. In particular, the oscillation frequency, the collapse and revival times are reproduced 
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astonishingly we... We flna..y discuss the corresponding ana.ytical expressions. It should 

be noted that the experimentally observed oscillation frequency in [1] of about 40ms follows 
directly from our formula. 



II. SEMICLASSICAL DESCRIPTION 



We will follow mainly Braun 12| and his discrete WKB method, as already applied to 
the double- well problem by Korsch et al. 13j. The two-mode Bose-Hubbard Hamiltonian 
can be written in the Schwinger spin representation by transforming to angular momentum 
operators Jx = ^ ia\d2 + d^ai ), Jy = (a|a2 — 0.20-1 ) and Jz 



^ ^otoi — 0002 ). With the 



ladder operators J+ = + iJy and J_ 



— ijy the Hamiltonian ([T]) becomes 



H = 2UJ: + 26 Jz -TlJ+ + l 



+ -UN^ - UN 
2 



(3) 



where N is the total particle number operator. For fixed a change from basis |n, N — n) 
to the angular momentum states |/, j) is useful, with I = ^ and j = . With Wj = 
2Uf-2Ul+2Uf+25i andpj = -T^yi{l + 1) - j(j - 1), the eigenvalues of the Hamiltonian 
are determined by an equation of the form 



PjCj_i + {wj - E)cj + Pj+iCj+i = 



(4) 



as discussed in 



12| . By introducing the "coordinate" operator 



'9] 



(note [IJ]), equation 



(jlj) can be written as a Schrodinger equation for the function Cj with eigenvalue E and 
Hamilton operator H = w{j) + p{j)e~'^^ + p{j + l)e*'^. In the classical limit the operators 
turn to canonically conjugate coordinate (p and momentum j (population imbalance), where 
turns out to be the phase difference between the two wells. Since p{j) is a slowly varying 
function of j in the classical limit (A^ — t- 00) one can replace both pj and pj+i by and 
one finds the Hamilton function 



H{j, 0) = w{j) + 2p{j + ^) COS. 



-UN^ -UN + 2Uf + 25] - 2T^{N/2Y 



j2 COS( 



(5) 
(6) 



which can ai. 
limit 
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so be found from the mean-field Gross-Pitaevskii functional in the two-mode 



16| . The classical dynamics of the population imbalance and the relative phase 
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(for (5 = 0) is then determined by Hamilton's equations of motion: 

dH 



dt 



d(t) 
d(t) dH 

dt d] ^WJW 



2Tv/(iV/2)2 -j2sin. 
2Tj cos 



3' 



(7) 



The rich dynamics in this "classical picture" have been studied by several groups 



ay 



18| , focusing on the differences between the classical and the quantum description of the 



dynamics 
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|-|2l| . Clearly a purely classical description cannot picture the collapses and the 
revivals of the population imbalance, but it is able to shed light on the transition from the 
tunneling to the self-trapping regime. Recently, the phase space region near the classical 



bifurcation was also investigated experimentally with ultracold Bosons 



22|. 



A. Semiclassical energy spectrum: Bohr-Sommerfeld quantization 

An analytical approach to the energ y sp ectrum relies on the WKB method following 
Brauu and others QQ. h Q. only the noninteract.ng case .s ...vestigated 
analyticalW. I. UM the authors eoncentrate on energies dose to the extremal points, 
and in [25[ the case of an attractive gas for the single value of A = 1 is studied. We here 
concentrate on the plasma oscillating regime and aim for solutions over the whole range of 
A < 1 to A > 1. 

For the Hamilton function (|6]) it is convenient to introduce two potential-energy curves 



= H{j, 7r) = TN + 2Uf + 2Ty/{N/2y-f (8) 



V-ij) = Hij, 0) = TN + 2Uf - 2Tv/(iV/2)2-j2 (9) 

such that the classically allowed energies lie in the region confined by the two potential 
curves and V~. The minimum energy is chosen to be V~{j = 0) = 0. The potential 
curves display the transition from the Rabi- to the Josephson regime very nicely, as shown 
in Fig. [H The energy eigenvalues change from a (non-harmonic) oscillator like spectrum 
for A < 1 to a spectrum with doublets for A > 1 due to tunneling, which can be seen 
from the potential curves. For A > 1, attains a local minimum which leads to doublets 
in the spectrum for energies E > V^{0). The deeper the minimum, the bigger this 
so-called Fock-fraction of the spectrum. Since we are interested in plasma oscillations, the 
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of 



FIG. 1: Potential curves V^, V and the energy eigenvalues £"„ for, from left to right, A = 0.1, 1, 
10 and TV = 30, T = 3. 

Fock-fraction will not be investigated here, but a semiclassical analysis along similar lines - 



if only more involved - is possible, see for instance 



26|. 



In the WKB approximation the eigenenergies En are obtained from the quantization 
condition 



S = SiE) = ^ <P{j)dj = 4 arccos ( ^^y^^^y— I = 2vr(n + -] 



(10) 



where n is the quantum number, and is determined by the Hamilton function (E]) at 
fixed energy E (recall that zero energy E = corresponds to S{E = 0) = 0). The integration 
limit j+ is the (positive) classical turning point as obtained from 



^{V+ij)-E)iE-V-{j)) = 0, 



(11) 



which leads to a quadratic equation in with solutions 

if)±iE) = ^ {{EU - {uj,l2f) ± ^{ujj2Y-EUujl/{2{l + K))) . (12) 



Recall that ujp = 2T\/1 + A is the plasma frequency. For the plasma oscillating regime 
the relevant turning point is j+. Note that j+ — )■ for — > 0, while (j^) approaches the 
negative constant (j^) — > —{ujp/{2U))'^ as E ^ 0. 

The integral in eqn. (fTOj) can be solved numerically and the results agree very well with the 
exact quantum results even for quite small numbers of particles as has already been noticed 



in [13[. It is impossible to solve the action integral analytically without approximation. As 

we aim at the plasma oscillation regime, we expand in powers of E. First, however, we take 

the derivative with respect to energy and rescale to find 

dS{E) ^ 2 f dX 
dE U\j4E)\J, ^(l-A2)(l + /t2(i^)A2) ■ ^ ^ 

with = Since Ofoi E ^ 0, and < A < 1, an expansion of (1+k^A^)~-^/^ 

in powers of k^A^ leads to a series in powers of E. The corresponding integrals d\ 
are known analytically. Finally, a systematic expansion of k^" and in E leads to 

dE Up a;3(l + A) + A)^ ' ^ ^ 

which is one of the important results of this paper. Apparently, the formal expansion in E 
is an expansion in the dimensionless parameter 

_UE _1 A ( E \ 

The expression on the right hand side clearly shows that our results are expected to be valid 
in the plasma oscillating regime E < V'^{0), irrespectively of the value of A. From a simple 
integration together with the Bohr-Sommerfeld-quantization condition (ITOl) we find 



^ ^ 2 Up u;|(l + A) w^(l + A)2 

In figure[2]we show examples of the spectrum for a wide range of values of A = 0.1, 1, 10, 
covering both the Rabi and the Josephson regime. Apparently, our approximation (fT6ll . 
including contributions up to third order in E, coincides with the numerically exact spectrum 
with high accuracy in the plasma oscillating regime {E < V^{0)) for all values of A. Clearly, 
the doublet structure in the Fock regime (high energy regime E > V^{0) in the right diagram 
of Fig. [21) cannot be captured by our series expansion (ITBl) . 



III. EXACT QUANTUM DYNAMICS OF THE POPULATION IMBALANCE 

To determine the tunneling dynamics, the Bose-Hubbard Hamiltonian ([1]) can be diago- 
nalized numerically for a finite number of Bosons. Using the eigenbasis {10™)}? the dynamics 
of \ipit)) is given by 

\m) = Yl cne-*^"*|0„) ; with c„ = (0„|^o) (17) 

n 
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FIG. 2: Comparison of the analytical (jl6p (red, solid line, including third order in E) and the 
numerically exact spectrum (blue squares) for, from left to right, A = 0.1, 1, 10 and for N = 30, 
T = 3. The vertical dashed lines illustrate the transition from the plasma oscillating regime to the 
self trapping regime at E = ^^'''(O) = 2NT. We see excellent agreement in the plasma oscillating 
regime. 



The time evolution of the population imbalance j = {fii — 712)/ 2 is then 



Kt) = (^(t)|j|^(t)) = 5^A_e-(^" 



-Em)t 



(18) 



with the matrix 



Anm = CnC^{(j)m\j\4>n) ■ 



(19) 



The dynamics of the population imbalance thus depends on the energy spectrum through 
the differences En — Em, and on the matrix Anm, which contains the initial condition and 
matrix elements {(j)rn\j\4'n) ■ 

Figure [3] shows the matrix Anm for increasing A, obtained from a numerically exact 
calculation. 
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FIG. 3: The matrix Anm for A = 0.1, 0.5, 1, 10 and for AT = 30, T = 3. 



Due to parity with respect to j, Anm is zero for an even number n — m, as can be seen 
in Fig. 131 = Ann = ^rm±2 = ^rm±4 + • • •• Amoug the non-zero matrix elements, there is a 



strong hierarchy, 

I » . . . , (20) 

in particular for small A, which will be important later. It is worth noting that in the limit 
t/ — i- (and therefore A — )■ 0, for fixed A^) the dynamics is well described by a harmonic 
oscillator. In that case it is easy to prove (the 0n(j) are Hermite polynomials) that only the 
Ann±i are in fact different from zero. 

Along the diagonals, the matrix elements Ann±k (with A; = 1, 3, 5, . . .) have a Gaussian- 
like n-dependence. This is due to the n-dependence of the overlap Cn = (f/'nlV'o), which will 
be discussed in the next section. By contrast, the n-dependence of the matrix elements 
{'Pn±k\j\'Pn) is weak. Thus, it is save to assume the form 

■Ann±k ~ ^nC^±kd'k 5 (21) 

with n- independent parameters dk ~ {4>n±k\j\4>n) (with the most relevant n), for which, 
following f l2U]) . we expect 

Mil > Msl > Msl > • • • (22) 



IV. SEMICLASSICAL DYNAMICS OF THE POPULATION IMBALANCE 

For a semiclassical evaluation of j{t) according to equations f|T8|) and f lT9|) we need semi- 
classical expressions for En — En±k and the overlap coefficients c„. While the spectrum was 
discussed in section |Tll we start here with the latter. 



A. Reflection principle 

The problem to find overlap integrals of an initial wavepacket ipoU) localized near j ^ jo 
with eigenstates \(j)E) of the Hamiltonian with potential V{j) = V^{j) is often encountered 



in molecular photo-dissociation |27|. The semiclassical solution (refiection principle) states 
that 

(^^\^o) = c■^Po(^^^^'^ (23) 

with some constant c. It is important to note that here the eigenstates are understood to 
be energy normalized, i.e. {(Pe\4>e') = S{E — E'), since in typical applications these are 
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scattering states. For the coefficients c„ we therefore find c„ = (^nlV'o) = y ^{4'e\'^o) ■ The 



normahzation condition 1 = ^„ |c„p ^ J dn|c„p yields c = 1/ \/V^JJq), and we get 



('^"I^O) ~ /777T-x V— ^0 T./.. N ■ (24) 



VnJo)^ dn V V'U, 

In our calculations, following the experiments, the initial wave function is prepared as the 
ground state of the tilted trap potential (achieved through the term 6{ni — 7x2) = 26j in the 
Bose- Hubbard Hamiltonian ([1])). In a harmonic approximation near the potential minimum 
of the tilted potential we find the Gaussian density 

\Mj)\' = ^e--^ (25) 



V2 



vro" 



with jo uniquely determined by the tilting strength S and 

^ , ^-P^°^^>^ (26) 
4 VI + A(l - (2jo/iV)=)-V' 

Clearly, the shape of the initial wave function determines the shape of the c„ as a function 

of n. On closer inspection of equation (l24l) . however, we observe that for the initial state 

( l25l) . due to the nonlinear relation between En and n, the coefficients c„ are gaussian in En 

but not in n. 



B. Population imbalance 

Having all the ingredients at hand we can now aim at a semiclassical expression for the 
dynamics of the population imbalance j{t) which we choose to write as 

jit) = ^ Ann-k exp{-i{En - En-k)t) + c.c. (27) 

n,k 

with k = 1,3,5,... taking into account the diagonal structure of Anm as discussed in the 
last section. Replacing Ann±k by expression (12T1) and using the Poisson summation formula 
we find 

oo 

j{t)= E E ^'mit) + c.c. (28) 



with 



(29) 
11 



This rather comphcated expression is readily simphfied by changing the integration variable 
from n to E. Further, as only very small (fc = 1, 3) are relevant (see equ. fl22|) ). it is 
safe to replace En — En-k ~ = 'i^W) neglect the fc-dependence in the reflection 
principle, i.e. c*_|_^ ~ c* . Finally, we replace 27rn = S{E) — n according to the semiclassical 
quantization rule (JTOl) . With t = kt we find 

2 

e"^e™^(-^). (30) 

This expression, together with equ. fl28|) is one of the main results of our paper. As we 
will see, even with further simplifications, the formula captures all essential details of the 
dynamics, allows for a thorough understanding of decay and revival dynamics, and, most 
importantly, is the starting point for analytical expressions. 



r 



dE 



^0 



E - VUo 
V'Uo) 



Due to the localization of the initial state ipo{j), the energy integration in fj30|) is 
confined to a relatively small interval near E ^ ^(jo)? which we assume to be in the 
plasma oscillating regime {E < ^+(0)). Therefore, for the evaluation of the overall phase 
mS{E) — 2TTkt/S'{E) we can rely on our semiclassical series expansions iHM and (fT6ll . 
With a Gaussian initial state as in f l25|) and expanding the overall phase up to second order 
around E ^ V{jo) allows us to take the Gaussian integral and leads us to the analytical 
result 

with T = kt. In the following we want to discuss the structure of this central result. The 
most important features are the plasma oscillations (a)p), their collapse (TcoUapse) and their 
revivals (Trev)- 

The phase if = (f{T,m) can be ignored for a qualitative discussion - it is a complicated 
expression and can be found in the appendix. Importantly, (f varies slowly with time and 
thus needs only be taken into account when quantitative agreement with exact calculations 
over extremely long time scales is sought. 

The parameter A = A{t, m) = t ■ — m ■ (expressions for the constants S,- and 
can be found in the appendix) describes an additional slow broadening and decay of the 
signal. As for the phase (p, the inclusion of A leads to quantitative agreement with exact 
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calculations as shown later, but need not be discussed further here. 

Thus we concentrate on the important plasma oscillations (cUp), their collapse (TcoUapse) and 
their revivals (Trev)- 

The analytical formula for the generalized plasma frequency for arbitrary A is 

Up = ujp{l - 2cig - ScsS-^) , (32) 

which is valid both in the Rabi and the Josephson regime. Here, ci = (1 + A/4)/(l + A) and 
C2 = (l + |^ + ^)/(l + A)^ are A-dependent numbers of the order of one and g = '^^2"'°^ is a 
dimensionless interaction parameter. We give a more elaborate discussion of this expression 
in section IVll 

For the revival time we find 

_ vr (1 + 2c,g) 

U{c, + 5c,gy ^^^^ 



and for the collapse time 



TcoUapse 2g AVo Up{c, + 5c,g) ' ^^^^ 



with 

AVo = crV'{jo)/Vijo) (35) 

being the width of the wavepacket in energy in units of the mean excited energy. Again, a 
more elaborate discussion of these results will be done in section I VII 



C. Simple Rabi limit 

In the well studied Rabi limit, i.e. when A ^ 1, our results simplify. In particular. 
Up — 7- Up, Trev — fj, and TcoUapse — > {2gAVoUp)~^. Moreover, only the main off-diagonal 
contribution k = 1 of the matrix Ann±k needs to be taken into account. Thus, in the Rabi 
limit, the dynamics of the population imbalance is governed by the simple expression 

J{t) = Jo M exp (-2uy{AVof (t - ^)'] , (36) 

a result that with an appropriate identification of the parameters can also be found in the 
literature 
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V. COMPARISON OF RESULTS 



Equation (13T|) describes the dynamics of the population imbalance without any free pa- 
rameter. The population imbalance oscillates with the generalized plasma frequency cUp, 
with roughly a Gaussian envelope of width TcoUapse (note that the parameter A contributes 
to the envelope, in particular for long times). The sum over m counts the revivals - the ini- 
tial collapse dynamics is captured by m = 0, the first revival corresponds to the contribution 
of m = 1, and so on. The sum over k takes into account further off-diagonal contributions 
in the matrix Ann±k which lead to small revivals (of the order of d^) at earlier times mTj-ev/k 
with A;— fold frequency. For A = 25, for instance, one can see tiny contributions of /c = 3 at 
one and two thirds of the full revival time in Fig. [71 

Figures HI O [6] and [7] show a comparison of the exact dynamics of the population im- 
balance, our analytical expression (13T]) (taking into account k = 1 only) and the simple 
expression for the Rabi limit f l36|) . for different values of A between 0.1 and 25. 
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FIG. 4: Comparison of the exact dynamics (top) of the population imbalance j with the improved 
semiclassical expression (j3ip (middle) and the expression for the Rabi limit (j36p (bottom) as a 
function of the dimensionless time t = u!pt/2iT for A = 0.1, T = 10, N = 100 and an initial jo = 20. 
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FIG. 5: Comparison of the exact dynamics (top) of the population imbalance j with the improved 
semiclassical expression (j3ip (middle) and the expression for the Rabi limit (j36p (bottom) as a 
function of the dimensionless time t = LOpt/lTi for A = 1, T = 10, N = 100 and an initial jq = 20. 

Obviously, our semiclassical expression ( I3T1) describes the exact dynamics almost perfectly 
over this huge range of values of A. By contrast, the simple expression (136|) is valid, indeed, 
for only very small values of A (A = 0.1), as expected. For increasing A the simple Rabi 
expression fails, as can be seen from Figs. [5] and |6l 




VI. DISCUSSION 

Having an analytical expression for the time evolution of the population imbalance allows 
us to discuss the dependence of the collapse- and revival time and the plasma oscillation 
frequency on the relevant parameters of the system. 
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FIG. 6: Comparison of the exact dynamics (top) of the population imbalance j with the improved 
semiclassical expression (j3ip (middle) and the expression for the Rabi limit (j36p (bottom) as a 
function of the dimensionless time t = Copt/lTi for A = 10, T = 10, N = 100 and an initial jo = 10. 

A. Plasma oscillation frequency 

The plasma oscillation frequency was found to be 

a)p = Wp(l - 2cig - 5025-^) (37) 

with ci = (1 + A/4)/(l + A), C2 = (1 + f + ^)/(l + A)2, and g = For very small A, 

the correct Up approaches the standard plasma frequency Up, since the constants Ci and C2 
tend to one in this limit, and the parameter g approaches zero: with V^(jo) ~ ^pilK'^^T) 
(harmonic approximation of the potential), it is worth writing the latter parameter in the 
form 

g ^ ^{2j,/Nf (38) 

which shows that g tends to zero linearly in A for fixed initial imbalance (Jq/N). However, 
for increasing A the correction terms in Up become more and more relevant, especially for 
large Jq, as can be seen from (|38|1 . Figure [9] shows a comparison of the plasma frequency 
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FIG. 7: Comparison of the exact dynamics (top) of the population imbalance j with the improved 
semiclassical expression (j3ip (bottom) as a function of the dimensionless time t = ojpt/2iT for 
A = 25, r = 10, iV = 100 and an initial jo = 10. 
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FIG. 8: Detailed comparison of the numerical exact initial collapse dynamics of j (green, dashed 
line) and the improved semiclassical analytical results (blue, solid line) for A = 25, N = 100, 
r = 10 and jo = 10 as a function of the dimensionless time t = ojpt/2-K. 



obtained from numerically exact results and the semiclassical expression ( 1321) as a function 
of A for different initial imbalance jo- H can be seen that the classical plasma frequency Up 
is only a good approximation for very small A, as expected. Especially for relatively large 
initial imbalance jo = 20 with = 100, the numerically exact plasma frequencies (circles) 
differ strongly from cUp, but are in very good agreement with the new semiclassical expression 
ojp (blue, solid line). Since A = 25 and an initial jo/N ^ 0.15 are typical experimental values 
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FIG. 9: Comparison of the numerically exact (symbols) and improved semiclassical analytical (j32p 
(lines) plasma oscillation frequency as a function of A for different jq: jo = 5: diamonds, dotted 
line, jo = 10: squares, dashed line, jo = 20: circles, solid line, for a total number of A = 100 
particles. The constant dash-dotted line indicates the simple plasma frequency ojp, valid only in 
the Rabi regime A <C 1. 

a. th,s discrepancy becomes by all means relevant. Sure enough, with the parameters gwen 
in [2|, our formula leads to 2tx jCj^ = 39ms - which is the experimentally observed value. By 
contrast, without our corrections one would find 271 /up = 30ms. 



B. Collapse time 



According to (ETj), the collapse time is given by 

= 2gAVoUp{c, + 5c,gy ^^^^ 

The expression in front of the brackets (which can be identified with the collapse time in the 
Rabi regime, i.e. for small A) can be approximated as N/{Aijjp{2jo/N))(T. Thus, assuming A 
and (jo/N) are kept fixed so that a is proportional to \/iV, the collapse time is proportional 



to \/N. This \/A^-behavior has been stated before in |25l. |29|. Our semiclassical formula 
shows, however, that this statement is only correct for the special case of fixed A and jo/N, or 
in the Rabi limit (A ^ 1). In all the other cases, the collapse time depends in a nontrivial 
way on through A and jo/N. Figures HI E], El [TJ and in detail Fig. M show that our 
semiclassical expression for the collapse time is remarkably reliable. 
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C. Revival time 



Following ( 13T|) . the revival time is 

T 

J. ro 



TT (1 + 2Ci( 



(40) 



U (ci + 5c2g) 

Most interestingly, in the Rabi limit it becomes independent of the number of particles 
and in fact independent of any other system parameter except the interaction strength U. 
The revival time was already discussed in 2l|] where it was found to be equal to An, with 
an interaction strength of | (considering the different definition of parameters) , which we 



confirm here, in the Rabi limit. Furthermore, it is stated in [21|, |25| that the revival time 
grows linearly with the number of particles A^. This is obviously true for those investigations 
with UN = const only, as can be seen from our expression (HOj) . However, note that even 
only slightly away from the Rabi limit, when A approaches or becomes greater than one, 
the constants Ci and C2 and the parameter g become relevant. This can be seen from Fig. 
[TOl Thus, for A > 1 no simple scaling law for the revival time exists. 




FIG. 10: Comparison of the numerically exact (symbols) and improved semiclassical analytical 
revival time (lines) as a function of A for different initial imbalance jq: jo = 5: diamonds, dotted 
line, jo = 10: squares, dashed line, jo = 20: circles, solid line, for a total number of = 100 
particles. The dashed-dotted line indicates the result ir/U of the Rabi limit, which is obviously 
only valid for very small A. 

The figure shows that for increasing A the exact revival times differ strongly from the 
revival time 71 /U predicted by the Rabi limit formula. On the other hand, it can be seen 
that the improved semiclassical expression (HOj) reproduces the exact revival times very nicely 
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even for A > 1. For increasing values of A the self-trapping fraction of the phase space is 
increasing as well, such that for large initial excitations, e.g. jo = 20 for = 100, the 
semiclassical analysis ceases to give reliable results. 

VII. CONCLUSION 

We applied semiclassical methods to the well-known two-mode Bose-Hubbard model, 
in order to investigate in detail BEC tunneling in a double-well trap. Within the plasma 
oscillation regime we found analytical expressions for the energy spectrum and the initial 
state agreeing nicely with numerically exact results. Employing the reflection principle and 
the Poisson summation formula led us to an analytical expression for the time evolution of 
the population imbalance of the Bose gas in the double well. This allows us to discuss the 
dependence of characteristic quantities of the dynamics, like plasma oscillation frequency, 
collapse and revival times, on the relevant system parameters. Remarkably enough, despite 
a wealth of publications on the two-mode model, such detailed understanding has not been 
achieved before. Finally, our generalized formula for the plasma oscillation frequency agrees 
perfectly well with experimental findings. Challenging as it may be, we hope that our 
predictions for collapse and revival times will be confirmed experimentally, too. 

Semiclassical methods are well suited to study the non-equilibrium dynamics of a Bosonic 
interacting many-body quantum system. For systems with more degrees of freedom, an 
explicitly time dependent approach might prove useful. 
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Appendix A: Parameters 

In order to complete the discussion of our semiclassical analytical result for the time 
evolution of he population imbalance (13T!) . we present here the definition of the remaining 
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parameters. The phase of the oscillation reads 



1 . A ( T — mTr, ^ ^ 



= -rrnp^ + - arctan A - — — , (Al) 

Z Z[i + A j \ -/collapse / 

where the dominantly m-dependent part is defined separately as 

cp^ = 2nV{l + Cig-l/i2V)) . (A2) 
V is the mean excited energy in units of the plasma frequency 

V = V{jo)/up . (A3) 

Furthermore, the quantity 

A = rE^ - m'Ern (A4) 

contributes to an overall slow spread and decay of the signal. It can be separated in a r- 
and a m dependent contribution with 

= lOc^iAVofg^Up (A5) 

and 

S^ = 47rci(A\/o)W- (A6) 

For completeness, we repeat the expressions for ci = (1 + A/4)/(l + A), C2 = 
(1 + f + f )/(l + A)2, g = and A^o = (tV'Uo)/V{jo) from section IVBl 
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